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Abstract 

In this paper we present a rigorous cost and error analysis of a multilevel estimator based 
on randomly shifted Quasi-Monte Carlo (QMC) lattice rules for lognormal diffusion prob¬ 
lems. These problems are motivated by uncertainty quantification problems in subsurface 
flow. We extend the convergence analysis in [Graham et ah, Numer. Math. 2014] to mul¬ 
tilevel Quasi-Monte Carlo finite element discretizations and give a constructive proof of the 
dimension-independent convergence of the QMC rules. More precisely, we provide suitable 
parameters for the construction of such rules that yield the required variance reduction for 
the multilevel scheme to achieve an e-error with a cost of 0{e~^) with 0 < 2, and in practice 
even 0 ft! 1, for sufficiently fast decaying covariance kernels of the underlying Gaussian random 
field inputs. This confirms that the computational gains due to the application of multilevel 
sampling methods and the gains due to the application of QMC methods, both demonstrated 
in earlier works for the same model problem, are complementary. A series of numerical exper¬ 
iments confirms these gains. The results show that in practice the multilevel QMC method 
consistently outperforms both the multilevel MC method and the single-level variants even 
for non-smooth problems. 


1 Introduction 


This paper gives a rigorous error analysis, together with numerical experiments, for a multilevel 
Quasi-Monte Carlo scheme applied to linear functionals of the solution of a typical model elliptic 
problem of steady-state flow in random porous media. This problem is of central importance in 
the development of efficient uncertainty quantification tools for subsurface flow problems. The 
random elliptic partial differential equation (PDE) reads 

— V • (a(x, t<;)Vri(x, a;)) = f{x), for x G D, w G Q , (1-1) 

where D is a bounded domain in for d = 1, 2 or 3, and Q is the sample space of a probability 
space (Q,.4, P), with cr-algebra A and probability measure P. A key feature is the coefficient 
a(-,a;), which is a lognormal random field on the domain D. 

In the context of flow through a porous medium, u is the hydrostatic pressure, a is the 
permeability and q := —aVu is the Darcy flux. This empirical relation between pressure and flux 
is known as Darcy’s law. When complemented by the conservation condition V • q = /, where 
f{x) is a deterministic source term, this leads to (|l.ip . 

In this paper, the uncertain permeability is assumed to take the form 


a(x, w) = a* (x)-b ao(x) exp 


j=i 


with Yj ~ A^(0,1) i.i.d. 


( 1 . 2 ) 


where a* and oq are given deterministic functions on D, satisfying o*(x) > 0 and ao(x) > 0. 
The sequence {tij} of nonnegative values is assumed to be enumerated in nonincreasing order. 
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accumulating only at zero, and the sequence {^j} is L 2 (-D)-orthonormal. If they correspond to the 
eigenvalues and eigenfunctions of the covariance operator of a correlated Gaussian random field, 
then the infinite sum under the bracket in (II.2p is known as the Karhunen-Loeve (KL) expansion 
of this Gaussian random field (see e.g. m)- 

For simplicity, we only study this problem subject to deterministic boundary conditions. 
In general, we may have mixed Dirichlet/Neumann conditions. Let the boundary F = dD be 
partitioned into two open, disjoint parts Tx> and Fj\^, and let n{x) denote the exterior unit normal 
vector to at X E F_/\^. Then we set 

u{x,-) = (t>T>{x) for X E Fx), (1.3) 

n(x) • (a(x, •)Vu(x, •)) = (t)j\f{x) for x E F_y\^. (1-4) 

For d = 2, 3, we assume D to be Lipschitz polygonal/polyhedral and each of Tx> and F^ to consist 
of the union of a finite number of edges/faces. 

Our goal is to obtain statistical information on certain linear functionals Q of the solution u to 
m; we write F := Q{u). In particular, we are interested in the expected value K[F] = E[^(u)] 
(with respect to the probability measure P). We need to perform several discretisation/truncation 
steps to obtain computable approximations to E[F]: 

(a) For a sample w, we employ a standard Galerkin finite element (FE) method with continuous, 
piecewise linear elements to discretise the solution to the PDE (11.11) on a family of simplicial 
meshes Fh parametrised by their mesh size h. We approximate entries of the element stiffness 
matrices by a one-point Gauss rule, that is, we evaluate the coefficient at the mid point of 
each mesh element. We denote the FE approximation on T^hy u^. 

(b) We truncate the KL expansion of log(a — a*) in (11.21) after a finite number of s terms; we 
denote the s-term truncated diffusion coefficient by and the corresponding PDE solution 
by Us- The FE approximation to (jl.ip on Th with a replaced by then reduces to a 
function Uh^g of s i.i.d. standard Gaussian random variables Y)-, j = 1,..., s. Denoting the 
approximation of F by Fh^s '■= G{uh,s), the expected value E[F] is then approximated by 

^ (i-s) 

where 4>{y) denotes the standard Gaussian probability density function. In porous media 
flow applications, the truncation dimension s is often very large. 

(c) The s-dimensional Gaussian integral in (II.5j) is then approximated by an iV-point quadrature 
rule, for example a Monte Garlo, sparse grid or Quasi-Monte Garlo rule, or by a multilevel 
variant (see below). 

In this paper the quadrature rules are derived from suitable Quasi-Monte Garlo (QMG) rules 
(i.e. equal weight rules on the s-dimensional unit cube), as we explain in the next section. The 
single-level variants of these rules, as estimators for (|1.5p . were analysed for the same model 
problem in the paper m (see also the earlier paper |23] for the uniform case). Much emphasis 
was placed there on the design of QMG rules that achieve dimension-independent error bounds 
with good convergence rates and under weak assumptions. 

Multilevel methods were introduced by [THl m- In the present context multilevel Monte 
Carlo (MLMC) estimators for (jl.5p (multilevel methods based on Monte Garlo integration) have 
attracted attention because of their capacity to reduce the cost without loss of accuracy. The 
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idea of using such multilevel estimators for the approximation of E[F] was established in [21 
and, for the lognormal case, analysed subsequently in [5l I34j . 

The multilevel method is based on a sequence of L + 1 FE approximations of increasing 
accuracy as ^ runs from 0 to L, with mesh diameters satisfying ho > hi > ... > h^. At level 
i we also truncate the KL expansion after S£ terms, with sq < si < • • • < sl- With the level 
£ approximation of our output functional F denoted by := we can write Fi as the 

telescoping sum 

L 

Fl = Fq + ^^{Fi — F>_i). (1.6) 

e=i 

Then by linearity of the expectation operator we have 

L 

EIFl] = E[To] + ^E[F, - Fe-i]. (1.7) 

1=1 

In the MLMC scheme each term is approximated by an independent Monte Carlo calculation, 
with a resulting gain in efficiency arising from the fact that the differences Fi — Fi-i on the higher 
levels, although more expensive to compute, have smaller variance and so require fewer Monte 
Carlo samples. 

In this paper, each of the L + 1 terms in (jl.7|) is instead approximated by a different QMC rule, 
where the number of quadrature points can again be chosen to decrease with £. For sufficiently 
smooth integrands, QMC quadrature rules offer the prospect of a higher accuracy for the same 
computational cost compared to standard Monte Carlo quadrature, or a lower cost for the same 
accuracy. Hence, the goal of this paper is to explore the combination of multilevel estimators 
and QMC methods by constructing and analysing a multilevel Quasi-Monte Carlo (MLQMC) 
estimator for the approximation of (fT3D . It was first observed in the context of stochastic 
differential equations in m that the two gains can be complementary. 

In the context of (ll.ip . single- and multi-level QMC FE approximations were analysed also in 
the recent papers PHIS], but for the simpler case of uniform and affine parameter dependence: 
in those papers the random variables Yj appeared linearly in the differential operator, and their 
values were assumed to be uniformly distributed on a bounded interval. The lognormal case 
considered here is technically more involved and the error bounds for the QMC rules developed 
here differ essentially from those for the uniform case. They require, for example, so-called 
“mixed regularity” of the solution of (II.5p . As shown here, this mandates stronger assumptions 
on the data than those required for MLMC or single-level QMC. The importance of this mixed 
regularity has already been recognised in |16] . In the present paper, we establish for the first time 
s-independent quadrature error bounds for MLQMC estimators and present detailed numerical 
experiments indicating that MLQMC methods can outperform single-level QMC and MLMC 
methods in terms of accuracy versus computational cost. Some numerical experiments have also 
been reported in [30]. 

The structure of this paper is as follows. Section [2| explains the mechanics of QMC methods, 
without entering into the question of approximation quality. Section [3| introduces the multilevel 
QMC method (MLQMC), establishes an abstract convergence theorem, compares the complexity 
of MLQMC to other estimators, and discusses practical aspects and a practical implementation. 
Section [D presents numerical results which confirm the theoretical results. All technical parts 
related to the necessary QMC convergence and construction theory are relegated to Section [5] 
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2 Quasi-Monte Carlo Quadrature 


Quasi-Monte Carlo quadrature rules are equal weight quadrature rules for integrals over the s- 
dimensional unit cube [0,1]^. For this reason we introduce a change of variables y = $7^(0) 
where := • • •, denotes the inverse cumulative normal distri¬ 

bution applied to each component of C £ [0,1]^. We then obtain from (jl.5h the expression 

= / F,,,($7nC))dC. (2.1) 

For the approximation of E[FX s] in a single-level scheme, we employ a specific kind of QMC 
quadrature rule, namely, the shifted rank-1 lattice rule given by 

Qs,N{Fh,s;A) := l^F,,,($7i(frac(^ + A))) , i = l,...,N, (2.2) 

where 2 G is the associated generating vector and A G [0,1]^ is the shift. The symbol 

frac(-) denotes the fractional part function, which is to be applied to every component of the 
s-dimensional input vector. For the general theory and fast construction of QMC lattice rules 
for the s-dimensional cube, see e.g., |20] as well as [291171 do]. For the particular case of integrals 
defined initially over see e.g., [2^ [28], 

The purely deterministic estimator ()2.2I) for E[F)j g] is biased. To remove this statistical bias 
we construct the associated randomly shifted lattice rule where the random shift A is uniformly 
distributed over [0,1]^. We then use the sample average of Qs,NiFh,s', over a hxed, finite 
number R of shift realizations as an estimator for E[F)i,,]. We arrive at 

1 ^ 

Qs,N,R{Fh,s) ■= J^'^Qs,N{Fh,s; ^k) , (2.3) 

^ k=l 

where Qs,NiFh,s', Afc) is defined in (|2.2p . k = 1,..., R. Now, let Ea[-] denote the expected value 
with respect to one or more random shifts. Since 

Ea[Q..„(F,,,;A)] = (trac()i+A))) dA 

1 ^ f 

= mT. T,,,($7i(A))dA = E[F,,,], 

the quantity in (|2.3p is an unbiased estimator for E[F)i^s]. However, (I2.3p is not an unbiased 
estimator for ElT], because the error arising from FE approximation and from truncation of the 
KL expansion of log(a — a*) cannot be removed by randomisation of (|2.2p . Specifically, the error 
analysis for randomly shifted lattice rules is carried out in terms of the root mean square error 
(RMSE) 

e{Qs,N,R{Fh,s)) := V^EA[(Q.,7v,fl(F;,,Q-E[F])'] . (2.4) 

Since the random diffusion coefficient a in (ll.ip is statistically independent of the random shift in 
the QMC quadrature rule, it is easy to see that in the single-level scheme we can split the RMSE 
as follows 

e{Qs,N,R{Fh,s))^ = ^A[{Qs,N,R{Fh,s) - + (E[F;,,, - F])\ (2.5) 
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The second term in (12.5|) is usually referred to as bias and can be decreased by choosing a 
fine enough FE mesh width h and by including a sufficiently large number s of terms in the 
KL expansion of log (a — a*), as discussed in [13]. The first term in (12.51) is the (shift-averaged) 
QMC quadrature error; it was analysed in detail in m where the crucial question of choosing 
the integer vector 2 in (|2.2I) was fully addressed. 


3 Multilevel Quasi-Monte Carlo Scheme 

Following the MLMC scheme, see mu and the subsequent MLQMC scheme for the uniform case, 
see |24] . we construct a multilevel Quasi-Monte Carlo estimator for M[F] by combining estimators 
of the form (|2.3I) on a hierarchy of levels. 

To define our multilevel method, let us assume that we have a nested sequence of FE spaces 
14 o, Vh ^, • • •, of increasing dimension and let Th ^, Th ^, ■ ■ ■, Thj^ be the corresponding sequence 
of shape-regular, conforming, simplicial meshes (i.e., simplicial partitions of the domain D for 
which intersections of any two d-simplices are are either empty, an entire side, or an entire face). 
We assume that the mesh diameters are strictly decreasing, i.e., . Furthermore, we 

include only the leading si terms in the KL expansion of logo on level subject to the condition 
si. < . The approximation of our output functional F that we obtain on level i is denoted 

by := Fh^^s^ as in (II.6p and for convenience we set E_i := 0. We can then write (II.7p as 


L 

e=o 

That is, the expected value of the output quantity of interest on the finest mesh is equal to the 
expectation on the coarsest mesh, plus a series of corrections, namely the expected value of the 
difference of quantities computed on consecutive FE meshes. We estimate the expected value 
K[Fi — Fi-i] on level i by means of the randomly shifted lattice rule estimator Qi := Qs^,Ni,Ri 
defined in ()2.3p and (12.2p . with quadrature points and Rg random shifts from a uniform 
distribution on [0,1]^^. The MLQMC estimator for E[F] then reads 



where := (frac [iZiN^ ^ + ^i,k)) and zi is the generating vector on level i (that will 

in general be different from level to level). 

Let us define the variance with respect to the shifts by 


Va[Q£(F> - F>_i)] = Ea[(QKF> - F>_i) - E[F> - F,_^]f] . 


Then, since each correction K[F£ — F£_i], £ = 0,..., L, is estimated using statistically independent 
random shifts, the RMSE of the MLQMC estimator satisfies 


L 

e(Q^L(^))2 _ EA[(Qr(E)-E[F])'] = A[QiiFe - F^.Q] + (E[Ei-F])^ (3.2) 

e=o 

The second term in (|3.2I) is the bias introduced by KL truncation and by FE approximation. It 
coincides with the second term of the single-level error in (12.5p for h = and s = s^. 


5 


3.1 Error versus cost analysis 

We now extend the cost analysis in [6l Thm. 1] to the MLQMC estimator defined in 

m- We aim at estimating the computational cost, denoted below by cost necessary 

to ensure that the RMSE in (j3.2p satisfied e (Q^^(E)) < s, as e 4 , 0. A similar extension of 
this abstract result has recently been proved in the context of multilevel stochastic collocation 
methods in [33]. However, our result here is tailored to MLQMC and includes the truncation 
error which was ignored in |33j . 

We assume the number of degrees of freedom := dim(V/i^), associated with the FE ap¬ 
proximation Fi := Ffif^se on level i = 0,..., L, satisfies 

Me ~ . (3.3) 

The assumption (|3.3I) includes quasi-uniform families of meshes and meshes with local refinement 
near corners or edges of the domain. 

Apart from the negligible post-processing cost to compute the quantity of interest, the cost 
of computing one sample on level £ is where denotes the cost 

of evaluating the s^-term truncation of the permeability field (11.21) at all quadrature points 
for each of the 0{hj'^) elements of the FE mesh, and denotes the cost of solving a sparse 

linear equation system with Me unknowns. We assume that 

~ ~ 7 >d. 

In the case of a robust (algebraic) multigrid solver, we have 'j = d +5, for arbitrarily small 5 > 0. 
In fact, the number of iterations for a robust multigrid solver typically grows only logarithmically 
with Me and the cost per iteration is 0{Me) (cf. |35] and the references therein). 

We will first state an abstract complexity theorem in which we make only very limited as¬ 
sumptions. To avoid having to treat the case £ = 0 separately, in the ensuing assumptions Ml - 
M3 we adopt the convention L_i := 1, s_i := 1, and recall that T_i := 0. 

Theorem 1 Suppose that lS.^[Qe{Fe — Fe-i)] = E(F£ — Fe-i) and that there are nonnegative 
constants a, a',/3,/3', 7 , A such that 

Ml. \E[Fl-F]\ < hl + sl^', 

M2. VA[Qt(T> - Fe-i)] < + (1 - , 

M3. cost(Q,(^t - Fe-i)) < Re Ne [seh^^ + hp), 

for all 0 < £ < L, and where 5.^. denotes the Kronecker delta. Then 

e{Qf^iF))^ < hr + + E {hh + (1 - and 

1=0 

L 

copQPiF)) < Y,ReNe{sehp + hp) . 

£=0 

Proof. The proof follows immediately from (|3.2p and the definition of cost(Q|^'^(E)). □ 

^Throughout the paper, the notation A < B indicates that there exists a constant c > 0 such that A < cB. The 
notation Ai^ B indicates that A < B and B < A. 
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We will now focus on a specific application of this theorem, with a fixed number of terms in 
the KL expansion. We assume that the sampling cost is the dominant part, which ultimately is 
the case with an optimal multigrid solver in the limit as the error tolerance goes to zero. We are 
not considering the case where the number of KL terms on the coarser levels is decreased, even 
though this may in some cases reduce the overall asymptotic cost of the multilevel algorithm, 
because it would lead to a very complicated complexity theorem and the analysis of Assumption 
M2 in Section [5] would become significantly more involved. 


Corollary 2 Letj < d+a/a' and let the assumptions of Theorem\^hold. If we choose ~ 2 
Ri = R and s£ = sl ^ some i? G N and for i = 0,..., L, then for any e > 0, there 

exists a choice of L and of Nq, Nl such that 


e{QT{F)f<e^ and cost{QfHF)) 


< 


g-2A-l/a' 

(log2 

^-2X-l/a'-{d-l3\)/a 


when fix > d, 
when PX = d, 
when PX < d. 


(3.4) 


Proof. Using the particular choices for h^, s^ and Ri and the assumption that 7 < d + a/a', we 
obtain 


+ and cost{Qf^{F)) < Y, Nehf^ . (3.5) 

£=o e=o 

Thus, a sufficient condition for the MSE to be bounded by a constant times is that each 
of the two terms in the above error bound is O(e^), which in particular leads to the choice 
2~^ ~ /iL ~ to bound the bias error, and thus 


L = 


1 


log2(e ^) + Cl 


a 


(3.6) 


for some constant ci G M that is independent of e. 

We now equate sampling and bias error to within a constant factor C 2 > 0, again independent 
of £ and of L To minimize the cost subject to this constraint, we consider the functional 


L / L 

g(iVo,...,iVi,/x) := AT, V'" + ^ - C2hi“ 

t =0 \ l =0 

where // is a Lagrange multiplier and where we treat Nq, ..., Nl as continuous variables. We look 
for its stationary point. This leads to the hrst-order, necessary optimality conditions 

=0 for £ = 0,..., L . (3.7) 

A 

L 

Y^I^'^4-a2hf =^- ( 3 - 8 ) 

£=0 


dg 

dNe 

dg 



Rearranging (13.7p . we see that NY~^^h^ i<^+h) -g independent of i. Therefore, the numbers of 
QMC points should be chosen according to 


Ne 


No 


f h^X 


for i = 1,..., L . 


(3.9) 
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A suitable choice for Nq can then be dednced from (|3.8p . Substituting (|3.9p into (13.8p and using 
the fact that h-o ~ 2^ = 1, we obtain ~ 2 ^“^ _ Since ^ 2~^, it follows 

from properties of geometric series that 


L L 

^(/3A-d)/(A+l) _ ^ 2^(rf_/3A)/(A+l) 

£=0 £=0 



^ El 


1 

L 

2L(rf-/3A)/(A+l) 


when /3A > d, 
when /3A = d, 
when /3A < d. 


(3.10) 


and hence 

No ~ 2^(2aA)^A_ 

Finally, we snbstitute (|3.9p and (|3.11l) into (13.51) and nse (|3.1UI) to bound that cost asymptotically, 
as L ^ 00 , by 


L 2^(2«'^+“/“') when/3A>d, 

COSt(Q^L(^)) < f^-a/a' ^^^f^iJ3X-d)/iX+l) < ^ ^E^aX+a/a') j^X +1 when/3A = d , 

e=o 2 ^( 2 “'’^+«/«'+<^-/ 5 '^) when / 3 A < d. 

The bonnd in (|3.4p then follows from (|3.6I) . i.e., using the relation 2^ ^ □ 

3.2 Discussion and comparison with other estimators 

First, let ns check the assumptions in Theorem [1] for the lognormal model problem (jl.ip . 

• We observe that Assnmption Ml relates only to the FE error and the KL truncation error, 
and is not specihc to MLQMC. It has been studied extensively in [U ISa [Ml EH]. The 
assnmptions on the data in Section [5l in particular on the regnlarity of the input random field 
a{-,uj) and of the functional Q, imply a = 2. For non-convex domains D, this reqnires special 
sequences of meshes and an analysis in weighted spaces (see Proposition 01 in Section 15.11 
which can also be used to bound the FE bias error). The valne for a' depends on the rate 
of decay of the KL eigenvalues. Under suitable regularity assumptions on the data, it was 
shown in 0] that, for Ganssian fields with Matern covariance and smoothness parameter ix 
(for a precise definition see Section 4), any a' < 2vjd can be chosen. 

• As shown in Section O the assumption that IEA[Qt(T£ — E£_i)] = E(F£ — E^-i) is satisfied 
for our randomised QMC rnles. 

• The main theoretical resnlt of this paper, postponed to Section |5l is to provide a proof of 

Assumption M2 for appropriate QMC rnles. We will see there that this assumption can 
usnally be satisfied for linear functionals, with (3 = 2a and with A G (1/2,1), for the case 
where sn = s^-i. The value of A, for a sufficiently good choice of the QMC rnles, depends 
on the parametric regnlarity of In particnlar, A can be chosen arbitrarily close to 

1/2 in the case of lognormal fields with Matern covariance and large enough smoothness 
parameter ix (as we discnss below). 

• Einally, if we use an optimal deterministic PDE solver, snch as mnltigrid, Assnmption M3 is 
also satisfied with 7 = (i + (5, for some 5 > 0, but typically S a/a' and thus 'j < d +a/a', 
as in Corollary [21 
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In practice, however, for the choices of parameters in Corollary [2] and assuming 7 ~ d, there 
is typically a critical tolerance e* > 0 such that for all e > e*. In that situation, we 

can drop the exponent — 1/a' in ()3.4I) for e > e^. Especially for d > 1, most practical choices for 
the tolerance e in applications lie above this critical tolerance e* > 0. We shall call the quantity 
obtained by dropping the —1/a' exponent the pre-asymptotic cost. Note however, that as seen 
in m, the QMC quadrature error also exhibits a pre-asymptotic behaviour. To obtain sharp 
bounds, the A in the pre-asymptotic cost should be replaced by the numerically observed effective 
rates l/Aefj < 1/A of the employed QMC rules. Note that the same is true for the single-level 
QMC estimator. There the cost is C)(^£- 2 A-i/a'-d/Q!j g^g ^ Q^^- 2 \ss-d/a^ ^ 

The analysis in mm of standard multilevel Monte Carlo (MLMC) methods for the lognormal 
case does not rely on the use of truncated KL-expansions. Isotropic input random fields a(-,w), 
such as those studied in Section 01 can be sampled in 0{h~‘^log{h~‘^)) operations via circulant 
embedding techniques (see, e.g., [H]). In that case, and so, with an optimal 

multigrid solver, the total cost on level i is 0{N£hj'^), for any 7 > d (for more details see 
Section 01 . Hence, assuming /3 7 ^ 7 , the cost of an optimal implementation of MLMC grows with 
^(g- 2 -max(o,( 7 -/ 3 )/a)) gg^£ ^ > d arbitrarily close to d. 

Nevertheless, for sufficiently large values of a' - typical for lognormal fields with Matern 
covariance and sufficiently large smoothness parameter 12 - we see that the presently proposed 
MLQMC estimator has significantly lower cost than, for example, MLMC estimators when A < 1. 
We will see in Section 0] that this holds in practice, even for values of the Matern parameter ly 
below the minimum required in the present convergence analysis. 


3.3 Practical aspects 

The formula (13.6h for L requires knowledge of the constant ci. When the error estimates are 
sharp, this can be computed a priori, as we do in our numerical experiments below. However, the 
FE discretization error, and thus the value of L, can also be estimated dynamically (i.e., without 
computing additional samples) from the estimates Qi{Fi — Fi_i), as for standard MLMC (see 

mm)- 

Like standard Monte Carlo estimators, randomised lattice rules also come with a simple 
variance estimator, namely the sample variance with respect to the random shifts, i.e.. 


Va[Q£(E,-F,_i)] 


1 

— 1 ) 


Rt 

[Qst,Ni{Fe. — E£_i; 


k=l 


Qse,Ni,Ri{Fi - Eg-i)]^ . (3.12) 


However, (on-the-fly) estimates for the rate of convergence 1/A of the lattice rule (or for its 
effective rate l/Aeg) are very unreliable, and thus the formulae (13.91) and (13.111) for the optimal 
values of Ni and Nq in the proof of Corollary [2] are of limited practical use. 

From a computational point of view, extensible lattice sequences or embedded lattice rules are 
useful, as they allow the results already calculated to be “recycled” when adaptively choosing the 
number of samples, see e.g., [laEiiiu]. To explore this “nestedness” property in practice, it is 
most convenient for the number of points Ni to be only powers of 2 (since then we always obtain 
complete lattice rules and do not need to be concerned about how the individual lattice points 
are ordered). A simple and effective algorithm that ensures this and does not require knowledge 
of A is presented in m- For completeness, let us recall the algorithm. To simplify notation, we 
define for £ = 0 ,..., L, := Va(Q£(T> - F>_i)) and Ci := cost{Qi{Fi - F>_i)). 


Algorithm 1 Let L = 0. 

1. Set Nl = 1 and estimate Vl using (j3.12p . 
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L 

2. While ^ V£ > e^, double on the level I for which the ratio V^/C^ is largest. 

£=0 

3. If the bias estimate is greater than e or L < 2, set L —>• L + 1 and go to Step 1. 


Note that this is a greedy algorithm that strives to equilibrate the profit, that is, the ratio of 
variance and cost, across levels. Thus, in the limit as e ^ 0, the numbers of samples Ng on the 
levels will be such that Vq/Cq ~ Vi/Ci Ri ■ ■ ■ Ki Vl/Cl- To show that this choice of Ni leads to 
the same overall cost for MLQMC as the theoretical algorithm in the proof of Corollary [2l let us 
assume that = vi (+ higher order terms), for some A > 0 and for some 0 < vi < that 

is independent of Ni. This is a stronger assumption than M2, but asymptotically it is satisfied 
for our QMC rules. Crucially, we do not require values of A, vi or j3 in the algorithm. 

We may also assume Ci = + lower order terms, where, at leading order, the “cost- 

per-sample” Kg ^ is independent of Ng. With these assumptions, we may set up a 

constrained optimisation problem, as in the proof of Corollary[21 minimising the total cost subject 
to the constraint in Step 2 of the algorithm on the total variance being less than However, 
here we write more abstractly 


g{No,, NL,fi) := '^Cg + fi i'^Vg - ^ 


e=o 


\£=0 


We ignore the higher and lower order terms in Vg and in Cg, respectively, treat the iVo ,... ,Nl as 
continuous variables again and differentiate g with respect to Ng and Jl to get 


dg 

dNg 

dfi 


Kg - ^Vg ^ 

= 0 . 
e=o 



Ng' = 0 


for i = 0,..., L , 


(3.13) 

(3.14) 


It follows from (j3.13p that Vg/Cg = \/fi, which is independent of i, and so the profit is indeed 
equilibrated across the levels for the optimal values of Ng. The fact that the asymptotic cost 
scales as in (13.41) can then be deduced as in the proof of Corollary [2l choosing 


Ng 


No 


\V0 Kg) 


A/(A+1) 


2 


where \x \2 '■= ^ 


that is, we round Ng up to the nearest power of 2. Substituting this into (13.141) . using (13.61) and 
the assumptions on vg and Kg, we can deduce that the expression for the optimal value for No is 
as in ()3.1ip (but rounded to the nearest power of 2). The bound on the cost follows as before. 

For standard multilevel Monte Carlo it is possible to compare this algorithm with the original 
algorithm in m that adaptively approximates the optimal choices of samples Ng, and we will 
see in Section 0] that Algorithm [1] achieves almost the same cost effectiveness as the original 
algorithm, even for fairly large e. 


4 Numerical results 

For all our numerical experiments we assume that the log-permeability loga(x,w) in (II.2p is a 
mean-zero Gaussian field with Mat&n covariance, that is, a* = 0, oq = 1 and {gj,^j) are the 
eigenpairs of the integral operator Pv{\x — x'\) v{x') dx' , with 
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where T is the gamma function and Ki, is the modified Bessel function of the second kind. The 
parameter is a smoothness parameter, is the variance and Ac is the correlation length scale. 
In practice, we will always truncate the sum in (ll.2p after a finite number of s terms. 

To compute the eigenpairs 1 < j < s, we discretize the integral operator above using 

the Nystrom method based on Gauss-Legendre quadrature on [0,1]'^ and then solve the resulting 
algebraic eigenvalue problem. 

The numerical results were obtained on a 2.4GHz Intel Core i7 processor in Matlab R2014b. 
4.1 Results in space dimension one 

We first consider problem (jl.ip in one dimension on D = (0,1) with homogeneous Dirichlet 
boundary conditions u{0,uj) = n(l,w) = 0 and source term / = 1. This problem is identical to 
the one studied in m Sect. 6]. For the discretization of the associated variational formulation 
on level i = 0,..., L we use piecewise linear, continuous FEs on a uniform simplicial mesh of 
width hi = ho2~^, where ho = for some £o £ such that Mi = _ i_ yyg generate 

samples of log a (and thus of a) at the midpoints of the intervals constituting the FE mesh using 
the KL expansion of log a with s terms, and approximate the entries of the stiffness matrix via 
the midpoint quadrature rule. The output quantity of interest is chosen to be F := u{l/3,u), 
i.e., the solution evaluated at x = 1/3. 

In order to have a nondimensional error measure for Q^^{F), our MLQMC estimator for 
E[E] with randomly shifted lattice rules, we define what is usually called the relative standard 
error in the statistical literature, that is 


erel{Qf^{F)) 


E[F] 


(4.2) 


We then study, for different tolerances e > 0, the computational cost to achieve a relative standard 
error erei(Q^'^(F)) < e and compare it to the cost to achieve the same relative standard error 
with standard MLMC, as well as with the single-level versions of both algorithms. In all the 
QMC estimators, for simplicity we use R = 16 random shifts and an embedded lattice rule with 
generating vector taken from the file [22l lattice-39102-1024-1048576.3600.txt]. (We remark 
that there is no theoretical justification to use this lattice rule for our problem here, however, 
numerical experiments from m indicated that such generic lattice rules do perform just as well 
as those specifically tuned to the problem.) 

We restrict ourselves to smoothness parameters ly > 1, where the numerically observed FE 
error is 0{h?) (independent of i^)Jl To estimate the bias error on the finest level L, we then 
assume the following upper bound (with uniform constants Cfe and Ctrunc): 

+ < Cfe+ Qrunc ^, (4.3) 

where F* is a reference solution computed with h* h and s* S> s (see [131 Sect. 2.4] for a 
justification). In Figure[I]we plot estimates of — F*]\ and of |E[F/i*^s — F*]\, for the case 

ofd=l,i/ = l, cr^ = l and for two different values of Ac. We also show bounds over the plotted 
range of h and s, for each of the two terms in (14.3p with the smallest possible values of Cfe and 
of Ctrunc- The expectations of these constants were estimated with 10^ MC samples and with 
h* = 1/1024 and s* = 500. We see that the rates of a = 2 = a' (for v = 1) in (14.3p are sharp. 

^Note that theoretically the FE error for point evaluations in one space dimension is 0{h? log |/i|) (cf. [32]), but 
we do not observe the log-factor in practice. 
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Figure 1: Estimates of the FE bias error — F*]\ (solid lines, left) and of the truncation 

error — F*]\ (solid lines, right), as well as the bounds in (14.3p for each case (dashed lines), 

for u = 1 and = 1. 


In our experiments, we then choose a particular sequence el ■= 2\/2C'fe^l> where L > 1 and 
CpE is the constant in (14.3p which we estimate as shown above for each problem. We choose a 
corresponding truncation dimension sl such that CtruncS^^’'^^ < Cfe^lj which implies 


SL ■ = 


Cbai 


with Cbal := (Ctrunc/CFE)'"/^'") , 


(4.4) 


and ensures that the total bound on the bias error in (j4.3p is less than / ^/2. We then run each 
of the estimators until the variance error is less than e^/2, thus ensuring a MSE (as defined in 
(j3.2p i of less than and a relative standard error (as defined in ()4.2p i of less than ei/|E[F]|. 

The numbers Ni of lattice points for the MLQMC estimator on each of the levels are chosen 
adaptively using the algorithm by Giles and Waterhouse [l2], given in Algorithm [1] in Section [3.31 
To estimate the variance Vi := Va(Q£(T£ — F£_i)) on each level, we use (I3.12p . As in Corollary[2l 
we choose si = on all coarser levels. For the cost on level £, we assume 

Cg := cost(Q£(F£ - Fi_i)) k. {2sl + 13) h~^ R . (4.5) 

This estimate is based on the fact (i) that the evaluation at the mid points of the mesh intervals 
of the coefficient in (11.21) . with a* = 0 , oq = 1 and with the sum truncated after sl terms, requires 
about = {2sL +1) operations; and (ii) that there are direct solvers for diagonally dom¬ 

inant tridiagonal systems (e.g., the Thomas algorithm) that achieve a complexity of 8 operations 
per unknown, leading to the cost estimate = 8{hJ^ + = 12hJ^. 

For the standard MLMC estimator we choose the same mesh and truncation parameters, h£ 
and Si, as for our new MLQMC estimator. The optimal numbers of samples are chosen 

according to the formula in the original paper m- This requires variance estimates for the 
differences Fi — Fi-i on each of the levels, which are obtained via the usual sample variance 
estimate with 10^ initial samples, updating the estimates as ^ oo on each level. The 

one-level variants are defined accordingly. 

In Figures [SHSl we plot the cost to achieve a relative standard error less than e with MLQMC 
and MLMC, as well as with the one-level variants QMC and MC, for = 1, v = 1,2, and 
Ac = 1,0.1. Red lines with circles correspond to the MC-based variants, while blue lines with 
diamonds correspond to the QMC-based estimators. The points on each graph correspond to the 
choices = 3 and L = 1,..., 4. The values of se are chosen according to (14.41) in each case. The 
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Figure 2: Cost to obtain a relative standard error less than e in the ID example. The covariance 
parameters are u = 2 and = 1, as well as Ac = 1.0 (left) and Ac = 0.1 (right), respectively. 
The estimates for Cbai are 0.76 (left) and 2.38 (right), respectively. 




Figure 3: Cost to obtain a relative standard error less than e in the ID example. The covariance 
parameters are u = 1 and ci^ = 1, as well as Ac = 1.0 (left) and Ac = 0.1 (right), respectively. 
The estimates for Cbai are 0.38 (left) and 1.78 (right), respectively. 


exception is the hardest test case {v = 1, Ac = 0.1), where we used L = 2,..., 5 and a variable 
number terms on level (. in MLQMC and in MLMC. The maximum number 

of KL terms included in that case is S5 = 456. In all test cases, we consistently see substantial 
gains for the MLQMC estimator, with respect to MLMC and QMC, even though the value of v 
is substantially smaller than our theory supports (see Remark fTOl ahead). 

For comparison, we show in Figures [2H3] also cost estimates for MLMC using circulant embed¬ 
ding which makes use of the Fast Fourier Transform (FFT) (magenta line, labelled ‘MLMC(FFT)’). 
Circulant embedding allows for efficient sampling at the quadrature points from isotropic random 
fields, such as the one studied here, without any truncation error (see e.g. M) and with a cost 
independent of sl. We assume that for the MLMC estimator with circulant embedding, the cost 
on level £ is 

^ ^ 2) h-^ iVf c < (68/9(£ + £ 0 ) - 248/27 + 12) {V2hi)-^ iVf ^ . (4.6) 

The factor y/2 in front of hg appears because there is no truncation error and thus the FE bias 
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Total Number of (Q)I\/IC Sample Points 



Figure 4: Left: MSE of the QMC/MC estimators for Fi — as functions of the total number 
of sample points, i.e. RN^ for QMC (solid blue lines) and for MC (dashed red lines), for 

V = = Xc = 1. Right: Total number of samples on each level to achieve a relative standard 

error less than e = 1.8 x 10““^ for the same example. 


error can be increased by a factor 2 to still achieve a MSE of e\ for the MLMC estimator. For the 
sampling of the coefficient we then assume the use of circulant embedding without padding [T3] 
- which doubles the number of unknowns in ID - and a split-radix FFT algorithm that requires 
^ n log2(n) — ^ n + C>( log2(re)) operations for vectors of length n [2T]. This is almost certainly 
underestimating the cost for circulant embedding, but, as we can see in Figures [2H3l the cost is 
still higher than that of our MLQMC estimator asymptotically. 

In Figure m we look at the particular case u = = Xc = 1, = 1/128 and sl = 49, and 

plot in the left figure the MSE of the QMC and the MC estimators for the expected values of 
the differences Fg — F/_i as the total number of sample points is increased (i.e., RN£ and 
respectively). We clearly see the faster rate of convergence with Ni ^ oo for the QMC estimators, 
which is almost optimal (i.e. the MSE is nearly 0{N^‘^)) even though = 1 is not sufficiently 
big for our theory in Section [5] to apply and even though in the construction of the QMC rules 
we did not use the weights derived there. We also clearly see the variance reduction from level 
to level (i.e., the offset between the lines), which does behave as theoretically shown in Section 5 
(i.e. roughly like 0{hj)). 

In Figure 0] (right) we plot for the same example the numbers of sample points on each of 
the levels. For MLQMC they were produced by Algorithm [H showing RN^, i.e. number of 
lattice points times number of shifts. For standard MLMC we show two sequences of numbers: 
those produced by the formula in the original MLMC paper labelled ‘MLMC(G)’, and those 
produced by Algorithm [T] with standard MC estimators on each level, labelled ‘MLMC(GW)’. 
We note that there are only very small differences in these final two sequences, conhrming our 
discussion in Section 1,8.,81 that Algorithm [D proposed in [12] can be used instead of the original 
algorithm to find the optimal sample distributions over the levels. The behaviour is the same for 
all other parameter values. 

4.2 Results in space dimension two 

We consider the problem (|l.ll) . (II. 2p with Mat&n covariance p,y in (|4.ip on D = (0,1)^ C 
At first we use again homogeneous Dirichlet conditions, i.e. F = Fx? and u(-,a;)|r = 0, and the 
source term / = 1. The output quantity of interest is the average of the solution u over the region 
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Figure 5: Left; Measured CPU times to calculate one sample of F) on level I for the 2D problem 
with Iq = 3 and s = 50. Right: MSE of the QMC/MC estimators for — F^-i], as functions 
of the number of sample points, i.e. RN^ for QMC (solid blue lines) and for MC (dashed 

red lines), for u = 1.5, = 1, Ac = 1, -^o = 3, L = 3 and sl = 27. 




F(w) := 


ID* 


Id* 


u{x,ui) dx. 


We discretise the associated variational formulation (spatially) using standard piecewise linear, 
continuous FEs on a sequence of triangular meshes obtained by taking a tensor product of each 
of the meshes in Section [Q with itself and by subdividing each of the squares of the resulting 
mesh into two triangles, thus leading to triangular elements of size := with 

ho := and Mi = ( 2 ^Mo _ (jggrggg gf freedom on level £ = 0,..., L. 

The finite element bias error and the truncation error are estimated as in ID. The choice of 
domain and functional guarantee that u(-,uj) G H^{D) (almost surely) and the FE and truncation 
errors converge as stated in (I4.3jl . for v > 1. Then, the number of KL terms is again chosen 
according to ()4.4p and si = sl on the coarser levels of the multilevel methods in all cases. For the 
average cost to compute one sample on each level, we use actual CPU-timings here (instead of 
FLOP counts). These were obtained using FreeFEM++ [T7| and the sparse direct solver UMFPACK 
[8]. The measured times to evaluate the KL expansion (with s terms) at the quadrature points 
(^perni) assemble and solve the sparse linear equation system (C|°^''*^) are shown in Figure 

[5] (left) together with the total time to compute one sample, for the case £o = 3, £ = A and s = 50. 
Finite element methods for dm) in two space dimensions allow, in the practical range of Mi 
considered here, for superior performance of sparse direct solvers as compared to, e.g., multigrid 
methods. Since we do not exploit the uniform grid structure in FreeFEM++ the cost in Figured] 
(left) is actually dominated by the FE system assembly, which scales like 0{hj‘^). We also note 
that <C for all our choices of below. 

In Figured] (right), we plot the MSE of the QMC and of the MC estimators for E,[Fi — Fi_i] 
as a function of the total number of sample points for the covariance parameters u = 1.5, = 1, 

Ac = 1, and for £o = 3, L = 3 and sl = 27. Again, we see the significantly faster and almost 
optimal convergence rate for the QMC estimators as iVg —>■ oo. 

In Figured] we plot again the cost to achieve a relative standard error less than e with all four 
estimators for two sets of covariance parameters. The points on each of the graphs correspond 
to the choices L = 1,..., 5 with = 3 (left) and L = 1,..., 4 with £o = 4 (right). We see 
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Figure 6: Cost to obtain a relative standard error less than e in the 2D example with homogeneous 
Dirichlet conditions, for v = 2.5, = 0.25, Ac = 1 and £o = 3 (left), as well as for ly = 1.5, 

cr^ = 1, Ac = 0.1 and £o = 4 (right). The estimates for Cbai are 0.55 (left) and 0.68 (right), 
leading to a maximum of S5 = 47 and S4 = 1106 KL terms on the finest mesh, respectively. 


similarly impressive gains with respect to MLMC and QMC in two dimensions, but we also see 
more clearly the influence of the smoothness parameter u. For the test case in the left figure, the 
numerically observed growth of the MLQMC cost is about 0(e“^'^^) over the range L = 1 to 4. 
For comparison, the costs for MLMC and QMC both show growths of 0(e“^'^) over the same 
range, while MC shows the expected 0{e~^) growth. 

As a final example, we consider the practically more interesting case of a 2D “flow cell”, 
that is, we solve the PDF (|l.ip in D = (0,1)^ with mixed Dirichlet-Neumann conditions. The 
horizontal boundaries are assumed to be impermeable, that is, (aVtt) • n = 0 for X 2 = 0 and 
X2 = 1. Along the vertical boundaries we specify Dirichlet boundary conditions and set u = 1, 
for xi = 0, and tt = 0, for xi = 1. We discretise this problem using the same sequence of meshes 
as above. Due to the Neumann conditions on the horizontal boundaries, the number of degrees 
of freedom in this problem is — 1 on level £ = 0,..., L. 

Here, the quantity of interest is the outflow through the right vertical boundary, i.e. 


F{u) 



a{x, uj) 


du{x, u) 
dxi 


dX2 . 

Xl=l 


As an approximation of this functional we use 


Fh^s{iyj) ■= - as{x,uj)Vuh,s{x,uj) ■ Vip{x) dx, 

Jd 

where ip denotes the FE function which is equal to one at all of the vertices of the right vertical 
boundary and is equal to zero at all other vertices (see [Ml section 3.4] for details). 

The numerical results for this problem are shown in Figure [71 In the left figure, we choose 
v = 2.5, = 1, Ac = 1. In the right hgure, we choose a set of parameters closer to the ones 

used in actual subsurface flow studies, namely u = I, = 3 and Ac = 0.3. In both cases £0 = 2. 
The points on the graphs correspond to the choices L = 1,..., 5 (left) and L = 2,..., 5 (right), 
respectively. The gains are again of the same order as above in both cases. In the smoother test 
case (left), the growth of the MLQMC cost is as low as between L = 3 and 5. 
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Figure 7; Cost to obtain a relative standard error less than e in the 2D flow cell example. The 
covariance parameters are v = 2.5, cr^ = 1, Ac = 1 (left) and v = 1, = 2>, Xc = 0.3 (right). 

The estimates for Cbai are 0.61 (left) and 0.0097 (right), respectively. The problem on the right 
is significantly more challenging. Please note the different range for e in the two figures. 

5 Mathematical analysis and construction of suitable QMC rules 

In the remainder of the paper, we present sufficient conditions on the data and on the FE spaces 
to verify Assumption M2 in the general MLQMC convergence result in Theorem [H as well as 
constructible QMC rules that achieve this. We will start in Section 15.11 bv addressing the spatial 
regularity and approximation orders for the FE function Uh^si'^u) in (|1.5p . making explicit the 
dependence on the parameter y in any constants that appear. Then we turn to the key estimates 
required for the MLQMC theory: bounds on the derivatives of the EE error with respect to the 
stochastic variables in certain weighted function spaces Wg which appear in the QMC convergence 
theory (see m and the references there) with constants that are independent of the truncation 
dimension s. These bounds correspond to “mixed derivative bounds”, appearing also in hyperbolic 
cross and other high-dimensional approximation methods [31311116], in that they require joint 
regularity of the random solution u{x,uj) with respect to the spatial as well as with respect to 
the stochastic argument. These estimates are proved in Section [521 and are used in Section [531 
to establish the MLQMC convergence rate estimates. 

5.1 Parametric formulation, spatial regularity and FE approximation 

As in [T3], we assume that, for d = 2, 3, D is a bounded, Lipschitz polygonal/polyhedral domain. 
Eor simplicity, we restrict ourselves to homogeneous Dirichlet data (p-z) = 0 and to deterministic 
Neumann data (/>_y G 77^/^(r_y) in dni). Then, the stochastic PDE (jl.l|l is (upon a measure-zero 
modification of the lognormal random field a in (|1.1I) L equivalent to the infinite-dimensional, 
parametric, deterministic PDE 

- V • (a(-,y) Vu(-,y)) = / in D, u|ri, = 0, n ■ a{-, y)V u{- , y)\r = (pjv, (5-1) 

with parametric, deterministic coefficient 

a{x,y) = a*(f)-bao(x)exp f ^ ^jix) yj j , (5.2) 

^ i>i ^ 
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where x £ D and where the parameter sequence y = {yj)j>i £ is distributed according 
to the product Gaussian measure JIq = A^(0,1). 

If y belongs to the set 

[/j, := E ^ bj |yj| < ooj C , (5.3) 

^ i>i ^ 

where the sequence b = {bj)j>i is defined by bj := ||CillL°°(_D) and is assumed to be in 

then the equivalence of (|5.1I) ~ (|5.2I) and (ll.ip - (ll.2p holds up to /I^j-measure zero modifications 
of the input random field. Due to the continuous dependence on the input random field, the 
parametric, deterministic coefficient a{-,y) in (15.21) and the parametric, deterministic solution 
u{-,y) of (15.ip will also differ only on a /Tg-nullset. If, moreover, b £ ^^(N), then the series (15.2p 
converges in L°°{D) for every y £ Ub, which we assume in what follows. 

To simplify the presentation, we assume a* = 0 and |rx)| > 0. We need the weak form of 
(j5.ip on the Hilbert space V = := {v £ H^{D) : u|ri, = 0}, with norm 

||u||y := ||Vu||i2p). 

As in m, we define for y £ Ub the parametric, deterministic bilinear form in V by 

£/{y;w,v) := / a{x,y)\/w{x) ■'S/v{x) dx, for all w,v£V. (5-4) 

J D 

We list properties of the parametric bilinear form £/{y;-,-) and of the weak solution u{-,y), as 
well as its FE approximation Uh{-,y), from [13]. For a proof see [131 Thm. 13]. 

Proposition 3 Assume that b £ £^(N). 

(a) The expressions 


a{y) := maxa(x, y) and 

a{y) := m.m.a{x,y) 

(5.5) 

x^D 

x^D 


are well-defined, JIq- measurable mappings from 

Ub to R which satisfy 


0 < a{y) < a{y) < oo 

for all y £ Ub ■ 

(5.6) 


(b) For every y £ Ub, the parametric bilinear form s^{y;-,-) : V x V ^ R defined in (15.4p is 
continuous and coercive in the following sense: 


£/{y;w,w) > d{y)\\w\\v 
£/{y]v,w) < a(y)||u||y||'u;||y 


for all w £ V, and (5-7) 

for all v,w £V . (5.8) 


(c) For every f £ L^{D) and every fix ^ o.nd with the (y-independent) linear 

functional 

Ffi{v) := / f{x)v{x)dx-\- / fix{x) v{x) ds , 

Jd Jtj^ 

the parametric, deterministic variational problem, to find u{-,y) £ V such that 

£/{y,u{-,y),v) = ^{v) for all v£V, (5.9) 


admits a unique solution u{-,y) £ V, for every y £ Ub- 
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(d) This parametric solution Ub 3 y u{-,y) G V is a strongly measurable mapping [with 
respect to a suitable a-algebra, cf. fT3^ ) which satisfies the bound 

lk(->y)llv ^ (||/IIl 2 (d) + ll</>Ar|li^i/ 2 (rAf)) /ora// yeUb, (5.10) 

[pointwise with respect to y). The implied constant depends on D, Tt> and Fw, but is inde¬ 
pendent of y. In particular, for any s G N, n(-, (y{i:s}; 0)) G V is well-defined^ measurable, 
and satisfies the above bounds uniformly with respect to s. 

(e) Restricting in to functions in the FE space 14 C V, there exists a unique, parametric 
FE solution Uh{-,y) G 14; for every y G Ub and 0 < h < 1, that also satisfies the bound 
()5.10l) . Replacing, in addition, the coefficient a{x,y) in (15.9p by the s-term truncated KL 
expansion a^{x,y), the corresponding s-parametric FE solutions Uh,s{'-, (y{i:s};0)) G 14 are 
uniquely defined for any yji.^} G M®, and satisfy, for (y{i:s};0) G Ub, the bound (j5.10ll 
uniformly with respect to h and to s. 

For the derivation of FE convergence rate bounds, we require additional spatial regularity: we 
assume in (15.2p 

a^ = 0, aoG and G . (5.11) 

With ()5.1ip . we may define the sequence 

:= 01“ max(||^j|4oo(^), |||V^j||4oc(o)), j = l,2,.... (5.12) 

Evidently, bj > bj so that <zUb C. We assume in what follows that 

b=(bfij>iGl\n) . (5.13) 

These conditions are satisfied under the provision of appropriate regularity of the covariance 

function of the Gaussian random field log(a — afi) in (11.111 (we refer to the discussion in m 
Rem. 4]). Also, for any b G f'^(N) C £^(N) we have ~ 

Proposition 4 Let us assume ()5.1ip and (|5.13p . and suppose we are given deterministic func¬ 
tions f G L‘^{D) and <j)jg' G Then the following results hold. 

(a) For every y G 17^, the series (15.211 converges in and 

a(-,y) g1T1’°°(Z1). (5.14) 


(b) The parametric solution map y i—>■ u{-,y) is strongly pLQ-measurable as a map from U-^ to 
the space 

W := {v gV ■. AuGL^{D)] , (5.15) 

and we have the a priori estimate 

l|An(-,y)| 42 p) < ri(y) {^Wfh^^D) + U^r\\H^/2(^^r)) > 


where 


Ti{y) := 


1 ||Va(-,y)||i,oo(D) 

d{y) a{yf 


for all yGU-^. 


(5.17) 


^As in for any finite subset u C N, we denote by (y^,; 0) the vector y G Ub with the constraint that yj — 0 
if j ^ u, and we use the shorthand notation {1 : s} for {1, 2,..., s}. 
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(c) There exists a sequence {V/i^}£>o of nested FE spaces of continuous, piecewise linear func¬ 
tions on conforming, simplicial meshes {Thi}£>o that satisfies the assumptions of Section\^ 
in particular, = dim(V/i^) ^ hf'^ and h^ ^ 2~^. The solutions Uhi{-,y) G defined in 
Proposition\^{e) satisfy the asymptotic error bound 

\\a^^‘^{-,y)'^{u-Uh,){-,y)\\LHD) ^ h(;T 2 {y)\\Au{-,y)\\L 2 (D), (5-18) 

where 

T 2 {y) ■.= d{yY/‘^ for all y^U-^. (5.19) 

The result holds verbatim for the FE solution Uh(,sei-, {y{i:s(}^^)) G 14^ of the si-term 
truncated problem. 

Proof, (a) This is a consequence of (|5.13p and V-G^^b) — 1 1311 Lem. 2.28]). 

(b) Since a{-,y) G n|ri, = 0 and (15.6h holds, for every y G the solution u{-,y) 

of (|5.ip also satisfies the following Poisson problem 

-Au{-,y) = f{-,y) := / ' [/ +Va(-,y) ■Vu{-,y)] inL^(D), (5.20) 

a{-,y) 

The bound p5.16p with Ti{y) defined in (I5.17P then follows from (15.101) . 

(c) The bound on ||An(-,y)|42(£)) in (b) together with the classical regularity theory for 
the Laplace operator on Lipschitz polygonal/polyhedral domains (see, e.g., [15]) implies weighted 
H‘^{D)-regularity of u{-,y) (with suitable weights near reentrant corners and edges) for non- 
convex D and full H‘^{D)-regularity for convex D. The existence of a sequence {14^}£>o that 
satisfies the assumptions of Section [3l together with 

inf lire —n||y < h£\\Aw\\]^ 2 (G>\ for all w , (5-21) 

then follows from classical FE results (for convergence bounds in weighted spaces see, e.g., |T]) 
and from the norm equivalence || At(;|42(£)) ~ ||w^||vK) for all w G W. The error bound in (|5.18p 
follows from an application of Cea’s Lemma m (in the energy norm) together with (|5.8p and 
(IOTP . □ 

Note that, for convex D and for homogeneous Dirichlet boundary conditions on all of dD, W = 
H^{D) and the family {Thg}i>Q can be constructed by uniform mesh refinement of an arbitrary 
conforming triangulation Fhq of D. 

5.2 Parametric Regularity 

As first observed in [24|, in the uniform case, the analysis of MLQMC methods for FE discreti¬ 
sations of PDFs requires estimates of the parametric solution map y i—?• u{-,y) in the regularity 
space W in ()5.15p . Here, we establish corresponding results for the lognormal parametric problem 
(j5.ip . (15.2p . In order to be able to draw upon our results in [T3|, we restrict the analysis to the 
particular case 

r^ = 0, r^ = r, V = H^{D), andV* = H-\D). (5.22) 

We denote by 5^ := {iz G : \u\ < oo}, where \u\ := (countable) set of all 

“finitely supported” multi-indices (i.e., sequences of nonnegative integers for which only finitely 
many entries are nonzero). For m,u G 5^, we write m < v '\i mj < vj for all j, we denote by 
u — m a multi-index with the elements Uj — mj, and we define (^) := nj>i (m ) • ^ sequence 

of non-negative real numbers (/3j)jgN we write (3‘' := nj>i ■ The following result is abstracted 
from the proof of [131 Thm. 14]. 
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Lemma 5 Given non-negative numbers {/3j)j^^, let o,nd be non-negative num¬ 

bers satisfying the inequality 




< 







for any ly G ^ {including v 


0 ). 


Then 




k<u 


^u—k 


where the sequence (A„)„>o is defined recursively by 


for all v G 


n—1 


Aq := 1 and An := ( . j Aj, for all n>l. 


i=0 


The result holds also when both inequalities are replaced by equalities. Moreover, we have 

nl 


An < 


(log 2y 


for all n > 0. 


(5.23) 


(5.24) 


Proof. We prove this result by induction. The case i/ = 0 holds trivially. Suppose that the 
result holds for all < n with some n > 1. Then for \u\ = n, we substitute m' = i/ — m in the 
recursion and use the induction hypothesis to write 


< 


E 

0^m'<u 


V 


V — m' 


r- K-rr^'V- 


< 


E 


jy 


u — m 




^2 ( jA|fc|+Bi^ 

k<u-m' ^ ^ 


Substituting k' = k -\- m', exchanging the order of summation, and regrouping the binomial 
coefficients, we obtain 


E (.:j E 

0^m'<u 


k' 

0^m'<k 


m'<k'<u' 


E (fc') E 


A\i,i_r„i\ + iu)i/ 


'I' 


k — m 


where 


E 

0 ^m<fc' 


k' 


m<k' 

m^k' 


k' 




E : V, . E E ; A. - E 


0 rn<k' 

\m\=i 


k' 




i=0 


k'\ 


Ai , 


which is equal to A|;;,/| as required. In the last step we used a simple counting identity (consider 
the number of ways to select i distinct balls from some baskets containing a total number of \k'\ 
distinct balls) 



(5.25) 


\m\=i 


The proof of ()5.24p then follows as in the proof of |13l Thm. 14]. 


□ 
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Theorem 6 For every f G L?‘{D), v' G ^ and y G with Ti{y) as in (15.171) . we have 


md''u{;y))\\L2^D) < \\f\\L^i^D)Tl{y)b^2\-\{\u\ + l)\. 


Proof. Throughout this proof, all estimates are for arbitrary y G t/5 C t/f, with the understanding 
that constants implied in ~ and < do not depend on y. For any multi-index i/" G 5^, we define 
(formally, at this stage) the expression 

9 u{-,y) ■■= V • (a(-,y) V((9'"'u(-,y))). 


Differentiation of order \u\ > 0 of the parametric, deterministic variational formulation (|5.9p 
with respect to y reveals that 

0 = £/{y,u{-,y),v) = / Vv{x) ■ d^'[a{x,y)Vu{x,y)) dx for all uGP. 

Jd 

The Leibniz rule d''{PQ) = Ylim<u (m) integration by parts imply 

V-d'^{aVu) = vY ^ =0 in P* . 


. .m, , 

m<u 

Separating out the m = u term yields the following identity in V' 

u 


V-(aV(a'^u)) = -V - f ^ iUd''~"^a)V{d^u) 


=g<^ 


E 

m<u 

E 

m<u 


m<u 

m^jy 

u 

m 


V ■ 


d‘'-^a 


{aV{d^u)) 


u \ f d‘ 


m 


' 9m d~ V 


d''- 


{aVid^u)) . 


In the last step we used the identity V ■ [pq) = pV ■ q + Vp ■ q. Due to (j5.6[) we may multiply 
g,^ by and obtain, for any \v\ >0, the recursive bound 


a g„ 


\\P{D) < X] 


m<u 

m^v 


V 

m 


L^{D) 


k ^^‘^9m\\L‘^{D) + 




V 


L°°{D) 


|aV2v(a-n)|k2p)^ . (5.26) 


By assumption, go = V ■ (a(-, y)Vn(-, y)) = —/ G L^{D), so that we obtain (by induction with 
respect to \u\ and using (15.lip and (I5.13p i from (I5.26P that y)yi/(-, y) G Lp‘{D), and hence 

from Proposition [31(a) that gv{-,y) G L^{D) for every u £ Thus, the above formal identities 
indeed hold in Lp‘{D). 

To obtain a bound on (|5.26p . we observe that it follows from the particular form of a in (|5.2I1 
that 

= {a — a^:)Y\{y/T'g for all v^m. (5.27) 

i>i 


22 



























Since we assumed a* = 0 in (15.lip , we then have 


d''-^a 


a 

L°°(D) 


i>i 


- n {v^Uj\\L<^{D) 

L^(D) j>i ^ 


= b" 


Moreover, using the product rule, we have 




k>l 


i>i 

j¥=k 


Due to the definition of bj in (15.121) . this implies, in a similar manner to (I5.28p . that 


V 


I I jy-m 

< \v — m\b 


L^{D) 


Substituting (|5.28p and (|5.29l) into (|5.26p . we conclude that 

II® ^ i ^ ^ ||a ^^‘^gm\\L^{D) 


m<u 

m^u 


where 


- E :) I" 

m<v' 




m<u 


I f TYl 1 -m fYt 

u — m\h A|^| b 


< A|,| b''^^ 
V^(y) \/b{y) 


In the first inequality in ^5.301) we used 

\\a^/^V{d^u)\\L^D) < A\m\b^ 


y* 




(5.28) 


(5.29) 


(5.30) 


(5.31) 


which was established in the proof of m Thm. 14]. In the second inequality in (I5.30h we used 
the bound bj < bj, for all j > 1, and the identity p5.25p to write, with n = ji^j. 


E 

m<u 


V 


m 


|i/-m|A|^l = 


n-l ^ ^ n-1 , K 

Y1 Y1 (^) (« - 0 A* = (J {n - i)^i =■■ A 

i=0 ^ ^ i=0 


i=0 m<i^ 
\m\=i 


Since Aq = ||a < Wfh^iD)/V^-iy), we now define 

« ^ T T*" 11 •^11^2 (D) ll/llv* 

Bjy := C'emb A|^| b —7=^ , where Cemb := sup 


■\/a(y) ’ feL2{D) II/IIl2(d) 

Then Aq < Bq and Bi, < for all v. We may now apply Lemma [5] to obtain 

-j- II/IIl2(^) 


A|fc| b^ Cemb A|^_fc| b 




(5.32) 
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Note the extra factor n — i m. the definition of A„ compared to in (|5.23l) so that A„ < A„. 
Using the bound in (j5.24l) . we have for all a < log 2 = 0.69 • • • , 


n—1 


n—1 




a 


n—i—\ 


n\ n\ „ n\ 

7-^ -TTT = - “ / TT^ - ae < — 

in — i — 1)! a” k\ 

^ ’ k=0 


n-1 


d 


2 = 0 ^ '' 2=0 

where the final step is valid provided that ae“ < 1. Thus it suffices to choose a < 0.567 • • • . For 
convenience we take a = 0.5 to bound (j5.32p . Using again the identity (I5.25p . we obtain 



kllliy — k\] = 




(|i^| + l)!. (5.33) 


Applying these estimates to (I5.32p gives 

l|a-'/VllL^(D) < Ce^b ”^^^ r2H(|^| + l)!. (5.34) 

va{y) 

Since • (aV(9'^u)) = a^/^A(cl‘^ri) + (Va • V(9'^n)), we have 

||a^^^A((9'^t6)||x,2(£)) < ||a ^^'^gu\\L'^{D) + II® (V®''^(^*^®))lli,2(£)), 

which yields 


V«(j/) II^(^'^®)IIl2(d) < 


la 


- 1/2 


'gv\\L'^{D) + 






a(f/) 


I® ' V((9 u)\\l2i^d) 


and in turn 


,1, / 11® ^ 9u\\lHd) , l|Va(-,y)||i»p) ||aV2v(a‘^u)||i2(^) 

11^(0' ®)IIl2(d) a -7==-^- z-r-. -7==-• 

\/aiy) a{y) \/a{y) 

Substituting (15.3411 and (|5.31l) into (I5.35p . and using A|j^| < 21*^1 |i/|!, as well as b'' < h'^, 
at 

II^(5‘'®)IIl2(D) < ^emh\\f\\L^{D)(^-Z^ + - -2l'^l (|i/| + 1)! . 


(5.35) 
we arrive 


This completes the proof. 


□ 


5.3 QMC convergence and design 


We first review the quasi-Monte Carlo theory that is essential for the QMC convergence rate 
estimates. We follow the setting and analysis in m Section 4] which, in turn, uses results from 
[28], see also the earlier references [36ll371[26l|25|. 

In our multilevel algorithm (|3.ip . for every level we apply a randomly shifted lattice rule 
Qi to the integrand — Fi_i which is multiplied by the product of univariate normal densities. 
Replacing F£—F£_i by a general function T" in s variables, we have the general integration problem 
J^s J~{y) 0^=1 c/)iyj)dy, with 4>iy) = exp(—y^/2)/\/^. The strategy in m is to consider a 
weighted function space with norm defined by 


ll-^ll 


w. 


(5.36) 



uC{l:s} 


aluljT 
-|u| dy^ 


(f/ui f/{l:s}\u) n ^iyi)<^y{l-.s}\u 


iG{l:s}\u 


nV'|(yi)d2/u, 

feu 


24 































where {1 : s} is shorthand notation for the set of indices {1, 2 ,..., s}, and denotes the mixed 
first derivative with respect to the “active” variables = {yj)j£u while ?/{i:s}\u = {yj)j&{i-.s}\u 
denotes the “inactive” variables. To ensure that the norm is finite for our particular integrand 
T = Fi — Fi-i, we follow [13] and choose the weight functions 

i’jiVj) = exp(-2aj |y|) , with 0 < Omin < aj < Omax < oo . (5.37) 


In Corollary[8]below, we will further impose the condition that aj > 9bj, with bj dehned by ()5.12p . 

A key ingredient in the analysis of lEI, see also [23l|2l], is to choose weight parameters 7u > 0, 
for every set u C N of hnite cardinality |u| < oo, such that the overall error bound does not grow 
with increasing dimension s. Such analysis makes use of the fact that the generating vector z for 
a randomly shifted lattice rule (see (I2.2p l can be constructed using a component-by-component 
algorithm to achieve a certain error bound, see [13p Thm. 15] or more generally |28] Thm. 8]. In 
particular, for F = Fi — Fi-i the result is that the variance (or the mean square error) of Qi is 
bounded by 


y^[Qi{Fi-F,_^)] < 



1/A 


btot(A'£)] 


-l/A 


\Fe — Fi 


£-1 


|2 

\Ws 


(5.38) 


feu 


for all A G (1/2,1], with 


/>,(A) := 2 


V^exp(a2/7*) 
7r2-2’?*(l — r/*)r/* 


C(A + ^), 


and 


V* ■= 


2A- 1 
4A 


(5.39) 


where := |{1 < z < — 1 : gcd(z, A^) = 1}| denotes the Euler totient function, and 

C(2^) := Yl^=i denotes the Riemann zeta function. Note that </?tot(A^) = A^ — 1 for A" prime 
and it has been verified that l/(/>tot(A) < 9/N for all N < 10^*^. Hence, from a practical point of 
view we can use 

V?tot(A) ^ N . (5.40) 

The best rate of convergence clearly comes from choosing A close to 1/2, but the advantage is 
offset by the fact that C(A -|- ^) —>■ oo as A —>■ 

To verify Assumption M2 in Theorem (TJ it remains to bound \\Fi — ™ (|5.38l) . Due 

to the triangle inequality. 


\\Fi '^sifWWai T ,S£_i) IIWs^ ) 

it follows from the next theorem and the subsequent corollary that M2 holds with /I = 4, in the 
case Si = S£_i. The remainder of the paper is then devoted to giving a choice of weights 7u that 
guarantees that the implied constant in M2 is independent of si. 


Theorem 7 Let s G N, /i > 0 and a* = 0, and let u be a general multi-index. Assume (15.lip 
and ()5.13p . Then, for every f G L'^{D) and for every Q G Lp‘{D)* with representer g G L?‘{D), 
we have for all y £ , 

\d‘'G{u{-,y) -Uhi-,y))\ < hf ||/||z, 2 p) 11^11^2(0) A(y)6''2l'"l(|iy| -^5)!, 
with an implied constant that is independent of h, f and g, as well as of y £ U-£, and with 

H{y) := T^{y)Ti{y) < oo, (5.41) 

where Ti{y) and T 2 {y) are as defined in (I5.17P and ()5.19p . respectively. 
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(5.42) 

(5.43) 


Proof. Let y E We define v^{-,y) E V and v^{-,y) E Vh via the adjoint problems 

£^{y,w,v^ {■,y)) = G{w) for all w&V, 

s^{y\Wh,v^{-,y)) = Q{wh) for all WhGVh. 

Due to Galerkin orthogonality for the original problem, i.e., 

■s^{y;u{-,y)-Uh{-,y),zh) = 0 for all Zh^Vh, 
on choosing the test function w = u{-,y) — Ufi{-,y) in ()5.42p . we obtain 

G{u{-,y) - Uh{-,y)) = £/{y,u{-,y) - Uhi-,y),v^{-,y) - v^{-,y)). 

From the Leibniz rule we have 

d''Q{u — Uh)= [ d''[aV{u — Uh) ■ V{v^ — v^)) dx 
Jd 

= j ^ {V{u-Uh)-V{v^-vl)) dx 

/ ^ {^Vd'^{u-Uh)-Vd^-'^{v^ -vl)dx 

^ m<u ^ k<m ^ 

- Uhi) ■ - uf)) dx. 


(5.44) 


k<m 


E 


u \ d'^-^a 


E 


m 

k 


m<i> ^ ' k<m 

Using the Cauchy-Schwarz inequality and (|5.28p . we obtain 


lirg(u-n,)i < E U 


m<u 


^ (^) k<m 


m 

k 


I a\'Vd'"{u — Uh)f dx j 


\ 1/2 


ID 


a\Vd^-'^{v^ - v^)\‘^ dx 


1/2 


< 


E 

m<i> 


V 


m 


■E 

k<m 


m 

k 


|ai/2va^(n - Uh)\\L2^n) - uf )||i2p) . (5.45) 


To continue, we need to obtain an estimate for — Uh)\\L2(^D)- Let X : U —)■ U 

denote the identity operator and let = 'Phiy) ■ ^ ^ Vh denote the parametric FE projection 
onto Vh which is defined, for arbitrary lu E U, by 

■^{y;'Ph{y)w-'w,Zh) = 0 for all ZhGVh- (5.46) 

In particular, we have Uh = VhU E 14 and 

Vl{y)=Vh{y) on Vh. (5.47) 

Moreover, since d^Uh E 14 for every fc E 5^, it follows from (|5.47l) that 

{T-Vh{y)){d'^Uh{y)) ^ d. (5.48) 

Thus 

||aL2v5^u - Uh)\\L2^D) = \W^^VVhd'^{u - Uh) + aL2v(X - Vh)d^{u - Uh)\\L2^D) 

< \\a^/^VVhd^{u-Uh)\\L2iD) + \\a^/^V{X-Vh)d'^u\\L2^^Dy (5.49) 
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Now, applying to (15.441) and separating out the £ = k term, we get for all Zh ^Vh, 


f aVd’^{u — Uh)-Vzhdx = — ^ ) f ^a)Vd^{u — Uh) ■ Vzhdx . (5.50) 


Choosing Zh = Vhd^{u — Uh) and using the definition (|5.46l) of Vh, the left-hand side of (I5.50p is 
equal to a |VP/j9^(u — u/i)p dx. Dividing and multiplying the right-hand side of (|5.5UI) by a, 
and using the Cauchy-Schwarz inequality, then leads to the bound 


[ a\V'Phd’^{u - Uh)\'^ dx 
J D 

dk-l^ 


E 


l<k 

t^k 


k 


a\Vd*^{u-Uh)\ dx] / a\Whd'^{u-Uh)\ dx 


’D 


L^{D) \JD 

Canceling one common factor from both sides and using (|5.28p . we arrive at 


a^/^VVhd'^{u-Uh)\\L2iD) < Y1 (5.51) 

£<k ^ ^ 

£^k 


Substituting (j5.5ip into ()5.49p . we then obtain 
\\a^/‘^Vd^{u - Uh)\\L‘2(D) 

f<k ^ ^ ^ V ^ ^ 

l^k A, Bfc 

Note that we have Aq = Bo- Now applying Lemma [5] with a = 0.5, Proposition SKc) and 
Theorem [6l we conclude that 


l|a'/W(n-n;,)||i2(^) < ^ A,,, ||aV2v(X - . 

£<k ^ '' 

< hT^iy) Y. (^)A|,|i-'||A(8‘-'u)||i.,o, 

£<k ^ ^ 

< h \\f\\L2i^D)Ti{y) Uy) Y, (fl 2l"l|€| 21''-^! {\k - £\ + 1)! 

£<k ^ ^ 

= h\\f\\L2^j,)n{y)T2{y)t2\^\ ( 1^1 + 2 )! ^ ( 5 _ 52 ) 

where Ti{y) and T 2 {y) are defined in (15.171) and (I5.19p . respectively, and where in the last step 
we used the identity 

j:Q)w!(ifc-<i+i)! = ffl^, 

which can be derived in the same way as (15.031) . 
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Since the bilinear form is symmetric and since the representer g for the linear 

functional Q{-) is in Lp‘{D), all the results in Section 15.11 hold verbatim also for the adjoint 
problem (j5.42l) and for its FE discretisation (|5.43p . Hence, as in (j5.52l) . we obtain 

(5.53) 

Substituting (I5.52p and (|5.53l) into (I5.45P yields 


WQ{U 


Uh)\ < 

E 


X 


m<u 


LHD) 

V 

m 




E 


k<m 


'm\-k (|fe| + 2)! -m-k ^\m-k\ (1^ 


k 


b 2i‘ 


k\ + 2)l 
2 


Using (I5.25|] we can obtain a similar identity to (j5.33p . 



{\k\ + 2)! (|m 
2 


fc| +2)! 
2 


(|m| + 5)! 
120 


Using again (j5.25p . with n = \iy\ we have 
^ \mj 120 

(^ + 5)! ^ ^ (i + l)(^ + 2)(z + 3)(i + 4)(i + 5)2- ^ (n + 5)! 

^\i) 120 120(n-U! “ 120 

2=0 2=0 


These, together with bj < bj for all j > 1, yield the required bound in the theorem. 


□ 


Now, to estimate the >V<j-norm of G{u—Uh), we need to bound its mixed first partial derivatives 
with respect to y = (yi,..., t/s, 0,0,...). The result in Theorem [7] was more general. In the 
following, we will only consider multi-indices u where each lyj < 1. As in the definition of the 
norm on Wg, we will use subsets u C {1 : s} of active indices instead of multi-indices. 


Corollary 8 Let clq := ao{x) and do '■= ao(x). For the weight functions ifj 

defined by (I5.37h with parameters aj > 9bj, we have 


||^(Ug Ufi^s 


\vvs 




E 

uC{l:g} 


[(H+5)!]- 

7u 


jGu 


with 


bj--= 


.t2 


exp(816j/2) ^{9bj 


aj - 9bj 


and 


K := \\f\\LfiD)\\9\\LfiD) + 


vrexp 



Proof. We begin by estimating H{y) defined in p5.41l) . It follows from p5.2p with a* = 0 that 

'Voo(x) 


Va{x,y) = a{x,y) 


ao{x) 


i>i 




leading to 


||Va(-,y)||ioo(£,) < d{y) 


||Vao||Loo(£)) 

ho 
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Since 1 + x < exp(x) for x > 0, we have 

||Vao||Loo(£)) 


H{y) < 1 + 


Oo 




i>i 


Kvf 

KvY 


<11+ 1 exp ( 2 ^ hj\yj \] ^ exp ( 7 ^ bj 


«o 


i>i 


i>i 


<11+ ^ Yl exp(% |yj|) . 


«o 


i>i 


Therefore, with := ||/||l 2 (£)) ||fi'||L2(D) (1 + il follows from Theorem [7] and 

the definition of the Wg-norm in (I5.36P that 




uC{l:g} 


[(|u|+5)!]2n,,u(45,^) 


7u 


(5.54) 


jG{l;s}\u 

leading to the univariate integrals 


n e^p(% \yj\)Y{yj) dy{i:g}\u I \yj\')Yj{yj) 

jeu 


1 < 


/ oo 

exp(% IvDYiy) dy 

-oo 


= 2 exp 


81 . j - 2 \ / j- n / 81.=-2 18 t 

< exp + ^6, 


and 


/ OO 

exp(18 6j |y|)V'|(y)dy 

-OO 


1 


aj - 9bj 


(5.55) 

□ 


These, together with (15.541) . then yield the estimate on the Wg-norm of ^(ti — Uh)- 

Theorem 9 For every f G L'^{D) and for every G G L?‘{D)* with representer g G L?‘{D), consider 
the multilevel QMC algorithm defined by (I3T]) with Si = s and Ri = R for all i = 0,..., L. 
Suppose that the sequence bj defined by (I5.12h satisfies 

ui 


b’^j < oo for some 0 < (? < 1 


i>i 

For each i = 1,..., L, let the generating vector for the randomly shifted lattice rule Qi be con¬ 
structed using a component-by-component algorithm fEEi, with weight parameters 

~ \ 2 /(l+A) 


7u : = 


(|u|+5)! 
120 


n 

j€u 


1 


\/QjW 


and a,- := - | 9bj + \l81bj + 1 “ ^ 


.t2 


1 


A := 


(5.56) 


in the weight functions (j5.37p . where bj is as defined in Corollary and 

2 ^ for some 5 G (0,1/2) when gG (0,2/3), 

2 ^ when gG(2/3,1). 

Let the generating vector for the randomly shifted lattice rule Qo constructed as in fT^ with A 
as defined in ()5.56l) . Then 

^A[Qi{Fi - Fi_i)] < D-,{X) R-^ hj , forall i = 0,..., L, 

where D-y{X) < oo is independent of s and £. 


(5.57) 
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Proof. First, let £ > 1. Using (|5.38|) and the triangle inequality, we obtain 


l/A 


yA[Qi{Fi — -p£-i)] < 


-11 


, 0^uC{l:s} j£u / 




-1/A 


Us '^/i£,s)llws “I" ||^('U'S llws^ • 


The bound in (|5.57|) now follows from Corollary [5] with 


D^{X) := 


|u|<oo jGu 


1/A 


|u|+5)!]- 

7u 



The fact that D-y{X) is finite can be verified following the same arguments as in the proofs of |13l 
Thm. 20 and Cor. 21]. 

Since by definition bj < bj and thus b G implies b G .f'?(N), the result for ^ = 0 follows 

from [T3l Thm. 20] with D-y{X) = C.y{X), defined in [T3l Eqn. (4.19)]. □ 


Together with (|5.40p , Theorem [9] shows that Assumption M2 of Theorem [T] holds with /3 = 4 
and A defined in (|5.56l) . 


Remark 10 As an example, let us consider the case where the KL expansion in ()1.2p arises from 
a Gaussian field with Matern covariance with smoothness parameter z^, as defined in Section HI 
We have from [131 Corollary 5] that fij < j-ii+’iu/d)^ Moreover, we see from the proof of [131 

Prop. 9] that ||V.^j||£,oo(£)) < for all d/2 + l<r<r<d + 2z/, allowing us to infer that 

bj < To ensure that the assumption < cc in Theorem [9] holds, we need 


q 




d/2 + l\ _ u-1 

d + 2v ) ^ d 


Therefore, a sufficient condition for the asumption to hold with g<lisz^>fi + l. A sufficient 
condition for q < 2/3 (and thus A = 1/(2 — 26)) is u > 3d/2 + 1. As we saw in Section (H these 
sufficient conditions do not seem to be necessary ones and we observe A ~ 1/2 even for much 
smaller values of u. 


Remark 11 Corollary [8] could be compared with [131 Thm. 16]. Unfortunately, there is a small, 
inconsequential error in m Eq. (4.17)]. The factors under the first product in [TH Eq. (4.17)] 
should be squared, and as a result, the denominator in m Eq. (4.11)] should also be squared. 
However, since this only amounts to the omission of a factor > 1 in the denominator, the estimate 
m Eq. (4.10)] is valid as stated. We have checked numerically that the weights [T3l Eq. (4.23)] 
with the adjusted formula for [13( Eq. (4.11)] lead, in all numerical experiments reported in [13j . 
to qualitatively the same results and therefore do not affect any of the conclusions drawn in [13| . 


Acknowledgments. Erances Kuo and Ian Sloan acknowledge the support of the Australian 
Research Council under the projects DP110100442, PT130100655, DP150101770. Robert Scheichl 
and Elisabeth Ullmann acknowledge support from the EPSRC under the project EP/H051503/1. 
Christoph Schwab acknowledges partial support by the European Research Council ERC and 
the Swiss National Science Eoundation SNSE during the preparation of this work through Grant 
ERC AdG247277. A large part of this work was performed during visits of Robert Scheichl and 


30 













Christoph Schwab to the School of Mathematics and Statistics, University of New South Wales, 

Australia, as well as the visit of Frances Kuo to the Department of Mathematical Sciences, 

University of Bath, UK. The authors thank Mahadevan Ganesh for suggesting an alternative 

proof strategy which led to quantitative improvements in the bounds of Theorem [6l and James 

Nichols for re-running all numerical experiments from [13] mentioned in Remark 111! 

References 

[1] T. Apel. Anisotropic Finite Elements: Local Estimates and Applications. Advances in Numerical 
Mathematics. Teubner, 1999. 

[2] A. Barth, C. Schwab, and N. Zollinger. Multi-level Monte Carlo finite element method for elliptic 
PDE’s with stochastic coefficients. Numer. Math., 119:123-161, 2011. 

[3] H.-J. Bungartz and M. Griebel. Sparse grids. Acta Numerica, 13:147-269, 2004. 

[4] J. Charrier and A. Debussche. Weak truncation error estimates for elliptic PDEs with lognormal 
coefficients. Stock. PDE: Anal. Comp., 1:63-93, 2013. 

[5] J. Charrier, R. Scheichl, and A.L. Teckentrup. Finite Element Error Analysis of Elliptic PDEs with 
Random Coefficients and its Application to Multilevel Monte Carlo Methods. SIAM J. Numer. Anal., 
51(l):322-352, 2013. 

[6] K.A. Cliffe, M.B. Giles, R. Scheichl, and A.L. Teckentrup. Multilevel Monte Carlo methods and 
applications to elliptic PDEs with random coefficients. Comput. Visual. Sci., 14:3-15, 2011. 

[7] R. Cools, F.Y. Kuo, and D. Nuyens. Constructing embedded lattice rules for multivariate integration. 
SIAM J. on Sci. Comput, 28:2162-2188, 2006. 

[8] T.A. Davis. Algorithm 832: UMFPACK V4.3—an unsymmetric-pattern multifrontal method. ACM 
Trans. Math. Software, 30(2):196-199, 2004. 

[9] J. Dick, F.Y. Kuo, Q.T. Le Gia, and Gh. Schwab. Multi-level higher order QMC Galerkin discretization 
for affine parametric operator equations. Preprint arXiv:1406.4432, Gornell University, 2014 (to appear 
in SIAM J. Numer. Anal. 2016). 

[10] J. Dick, F. Pillichshammer, and B.J. Waterhouse. The construction of good extensible rank-1 lattices. 
Math. Comp., 77:2345-2374, 2008. 

[11] M.B. Giles. Multilevel Monte Carlo path simulation. Oper. Res., 56(3):607-617, 2008. 

[12] M.B. Giles and B.J. Waterhouse. Multilevel quasi-Monte Garlo path simulation. Radon Series 
Comp. Appl. Math., 8:1-18, 2009. 

[13] LG. Graham, F.Y. Kuo, J.A. Nichols, R. Scheichl, Gh. Schwab, and I.H. Sloan. Quasi-Monte Carlo 
finite element methods for elliptic PDEs with lognormal random coefficients. Numer. Math., 131:329- 
368, 2015 

[14] LG. Graham, F.Y. Kuo, D. Nuyens, R. Scheichl, and I.H. Sloan. Quasi-Monte-Garlo methods for 
elliptic PDEs with random coefficients and applications. J. Comput. Phys., 230:3668-3694, 2011. 

[15] W. Hackbusch. Elliptic Differential Equations: Theory and Numerical Treatment, volume 18 of 
Springer Series in Computational Mathematics. Springer, 2010. 

[16] H. Harbrecht, M. Peters, and M. Siebenmorgen. Multilevel accelerated quadrature for PDEs with 
log-normal distributed random coefficient. Preprint 2013-18, Math. Institut, Universitat Basel, 2013 
(to appear in Math. Comp. 2016). 

[17] F. Hecht. New development in FreeFem-f-|-. J. Numer. Math., 20(3-4):251-265, 2012. 

[18] S. Heinrich. Multilevel monte carlo methods. In Lecture Notes in Large Scale Scientific Computing, 
number 2179, pages 58-67. Springer-Verlag, 2001. 


31 


[19] F.J. Hickernell and H. Niederreiter. The existence of good extensible rank-1 lattices. J. Complexity, 
19:286-300, 2003. 

[20] F.Y. Kuo J. Dick and l.H. Sloan. High-dimensional integration - the Quasi-Monte Carlo way. Acta 
Numerica, 22:133-288, 2013. 

[21] S.G. Johnson and M. Frigo. A modified split-radix FFT with fewer arithmetic operations. IEEE 
Trans. Signal Processing, 55(1):111-119, 2007. 

[22] F.Y. Kuo. Lattice rule generating vectors, web.maths.unsw.edu.au/~fkuo/lattice/index.html 

[23] F.Y. Kuo, Ch. Schwab, and l.H. Sloan. Quasi-Monte Carlo finite element methods for a class of elliptic 
partial differential equations with random coefficient. SIAM J. Numer. Anal., 6(50):3351-3374, 2012. 

[24] F.Y. Kuo, Ch. Schwab, and l.H. Sloan. Multi-level Quasi-Monte Carlo finite element methods for a 
class of elliptic partial differential equations with random coefficient. Found. Comput. Math., 15:411- 
449, 2015. 

[25] F.Y. Kuo, l.H. Sloan, G.W. Wasilkowski, and B.J. Waterhouse. Randomly shifted lattice rules with 
the optimal rate of convergence for unbounded integrands. J. Complexity, 26:135-160, 2010. 

[26] F.Y. Kuo, G.W. Wasilkowski, and B.J. Waterhouse. Randomly shifted lattice rules for unbounded 
integrals. J. Complexity, 22:630-651, 2006. 

[27] M. Loeve. Probability Theory, volume 11. Springer-Verlag, New York, 4th edition, 1978. 

[28] J.A. Nichols and F.Y. Kuo. Fast component-by-component construction of randomly shifted lattice 

rules achieving convergence for unbounded integrands in in weighted spaces with POD 

weights. J. Complexity, 30:444-468, 2014. 

[29] D. Nuyens and R. Cools. Fast algorithms for component-by-component construction of rank-1 lattice 
rules in shift-invariant reproducing kernel Hilbert spaces. Math. Comp., 75:903-920, 2006. 

[30] P. Robbe, D. Nuyens, and S. Vanderwalle. A practical multilevel quasi-Monte Carlo method for elliptic 
PDFs with random coefficients. In Master Thesis “Een parallelle multilevel Monte-Carlo-methode 
voor de simulatie van stochastische partiele differentiaalvergelijkingen” by P. Robbe, June 2015. 

[31] Ch. Schwab and C.J. Gittelson. Sparse tensor discretizations of high-dimensional parametric and 
stochastic PDEs. Acta Numerica, 20:291-467, 2011. 

[32] A.L. Teckentrup. Multilevel Monte Carlo Methods for highly heterogeneous media. In C. Laroque, 
J. Himmelspach, R. Pasupathy, O. Rose, and A.M. Uhrmacher, editors. Proceedings of the 2012 
Winter Simulation Conference. WSC, 2012. 

[33] A.L. Teckentrup, P. Jantsch, C.G. Webster, and M. Gunzburger. A multilevel stochastic colloca¬ 
tion method for partial differential equations with random input data. SIAM/ASA J. Uncertainty 
Quantification, 3:1046-1074, 2015. 

[34] A.L. Teckentrup, R. Scheichl, M.B. Giles, and E. Ulhnann. Further analysis of multilevel Monte Carlo 
methods for elliptic PDEs with random coefficients. Numer. Math., 125(3):569-600, 2013. 

[35] P.S. Vassilevski. Multilevel Block Factorization Preconditioners. Springer, 2008. 

[36] G.W. Wasilkowski and H. Wozniakowski. Complexity of weighted approximation over K^. J. Approx. 
Theory., 103:223-251, 2000. 

[37] G.W. Wasilkowski and H. Wozniakowski. Tractability of approximation and integration for weighted 
tensor product problems over unbounded domains. In K.T. Fang, F.J. Hickernell, and H. Niederreiter, 
editors, Monte Carlo and Quasi-Monte Carlo Methods 2000, pages 497-522, Berlin, 2002. Springer. 


32 


Prances Y. Kuo 

f.kuoOunsw.edu.au 

School of Mathematics and Statistics 

University of New South Wales 

Sydney NSW 2052 

Australia 


Robert Scheichl 
R.ScheichlObath.ac.uk 
Department of Mathematical Sciences 
University of Bath 
Bath BA2 7AY UK 


Christoph Schwab 

Christoph.schwabSsEun.math.ethz.ch 

Seminar fiir Angewandte Mathematik 

ETH Zurich 

Ramistrasse 101 

8092 Zurich 

Switzerland 


Ian H. Sloan 

i.sloanOunsw.edu.au 

School of Mathematics and Statistics 

University of New South Wales 

Sydney NSW 2052 

Australia 


Elisabeth Ullmann 
elisabeth.ullmannOma.turn.de 
Department of Mathematics 
Technische Universitat Miinchen 
BoltzmannstraBe 3 
85748 Garching 
Germany 


33 



